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Abstract — The  ability  to  precisely  determine  the  location  of 
radar  transmitters  can  be  crucially  important  in  maintaining 
domain  awareness.  This,  however,  may  be  problematic  with 
traditional  methods  when  used  with  a  distributed  network  of 
disparate  sensors.  A  novel  geolocation  technique  for  circularly 
scanning  radar  transmitters  is  introduced.  This  technique  uses 
the  differenced  central  times  of  arrival  (DCTOA)  of  the  main 
beam  as  an  observable.  The  solution  for  the  transmitter’s 
position  and  scan  rate  are  given  using  a  weighted  least  squares 
approach  as  well  as  a  particle  swarm  optimizer.  Experimental 
results  show  this  technique  is  able  to  locate  a  radar  transmitter 
within  11  meters,  while  maintaining  minimal  complexity.  This 
technique  has  the  advantage  of  requiring  orders  of  magnitude 
less  timing  synchronization  among  receivers,  an  order  of  mag¬ 
nitude  less  data  transfer,  and  it  does  not  require  simultaneous 
illumination  of  receivers. 
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1.  Introduction 

The  ability  to  accurately  locate  radar  sources  in  a  timely 
manner  is  crucial  in  maintaining  maritime  and  other  domain 
awareness.  Traditional  methods  of  precision  geolocation  re¬ 
quire  multiple  receivers  to  be  time  synchronized  on  the  order 
of  nanoseconds  and  often  require  simultaneous  illumination. 
Further,  receivers  must  exchange  data  that  often  occurs  at 
the  pulse  rate  of  the  radar.  This,  however,  may  present  a 
challenge  when  a  distributed  network  of  diverse  sensors  is 
used  that  have  limited  timing  and  data  exchange  abilities.  A 
method  of  precision  geolocation  with  less  stringent  timing 
and  data  exchange  requirements  would  enable  a  disparate 
network  of  sensors  to  accurately  estimate  a  radar  transmitter’s 
location,  as  well  as  expand  the  state  of  the  art  of  radar 
geolocation  techniques. 

There  are  a  number  of  methods  currently  used  for  radar 
geolocation.  Reference  [1]  provides  a  general  survey  of 
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techniques,  while  Reference  [2]  discusses  the  theory  of  these 
techniques  in  detail.  Many  of  these  techniques  require  the 
same  radar  signal  to  be  observed  by  multiple  receivers. 
A  common  example  of  this  is  Time  Difference  of  Arrival 
(TDOA),  where  the  difference  in  signal  time  of  arrival  is 
measured  between  two  spatially  distributed  receivers.  De¬ 
tails  of  this  technique  can  be  found  in  Reference  [3].  This 
technique  may  be  problematic  with  arbitrarily  distributed 
receivers  as  it  requires  simultaneous  illumination  and  precise 
timing  synchronization  between  receivers  (as  solution  error  is 
proportional  to  timing  error  times  the  speed  of  light).  Other 
common  techniques  have  similar  requirements  but  observe 
different  signal  characteristics.  For  example,  the  difference 
in  the  phase  of  the  signal  between  two  receivers  may  be  ob¬ 
served  in  interferometry  or  the  Difference  in  the  Frequency  of 
Arrival  (FDOA)  may  be  used  to  estimate  a  radar  transmitter’s 
location.  These  techniques  have  similar  timing  constraints 
and  will  further  require  data  to  be  shared  between  receivers  at 
the  pulse  rate  of  the  radar,  which  may  be  prohibitive. 

Additionally,  the  doppler  shift  of  the  signal  may  be  observed, 
if  there  is  relative  motion  between  the  transmitter  and  the 
receiver,  and  can  be  used  to  estimate  a  transmitter’s  location. 
This  technique  may  be  problematic  in  practice  as  it  gener¬ 
ally  requires  information  about  the  emitted  frequency,  which 
may  not  be  readily  available  and  may  vary  over  time  in  an 
unknown  fashion.  Reference  [4]  contains  further  discussion 
of  this  method. 

The  last  common  technique  is  based  on  the  Angle  of  Arrival 
(AOA)  of  the  signal.  Reference  [5]  as  well  as  Reference  [2] 
give  ample  discussion  of  this  method.  This  technique  uses 
multiple  measurements  of  the  AOA  of  a  signal  to  triangulate 
its  source.  This  method  does  not  require  precision  timing 
synchronization  of  receivers  or  high  rates  of  data  exchange. 
However,  as  determining  the  AOA  of  a  signal  to  sufficient 
accuracy  may  be  challenging,  precise  geolocation  results  are 
not  often  achievable, 

A  novel  method  has  been  developed  that  uses  differenced 
times  of  arrival  of  maximum  signal  strength  across  multiple 
receivers  to  locate  a  circularly  scanning  radar  transmitter. 
The  method  allows  for  passive,  precision  geolocation  without 
precise  timing  requirements  or  prohibitive  data  exchange 
volumes.  The  difference  between  signal  times  of  arrival 
is  related  to  the  angle  the  radar  has  swept  out  between 
observation  events,  thus  the  radar  transmitter  location  may 
be  calculated  to  be  along  a  line  of  constant  angle  between  the 
two  receivers.  Given  sufficient  data,  the  radar  transmitter’s 
location  may  be  estimated.  This  method  has  the  advantage  of 
being  independent  of  platform  and  tolerant  to  timing  error  on 
the  order  of  hundreds  of  microseconds. 

This  technique  has  been  developed  and  tested  at  the  US  Naval 
Research  Laboratory  (NRL),  and  has  been  developed  by  other 
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researchers  in  References  [6],  [7]  and  [8].  Major  work  at  NRL 
include  field  tests  in  2006  and  2009,  the  latter  of  which  is 
documented  in  Reference  [9].  NRL  has  been  issued  a  United 
States  patent  for  this  algorithm,  [10]. 

A  complete  description  of  the  methods  used  to  reliably  imple¬ 
ment  this  technique,  called  SearchLight,  are  presented.  The 
development  of  the  basic  observable  is  presented  along  with 
the  necessary  algorithm  to  allow  the  technique  to  be  used 
for  an  arbitrary  receiver  configuration.  Next,  methods  for 
estimating  the  radar  transmitter’s  position  are  presented  for 
both  the  linearized  and  the  nonlinear  system.  The  nonlinear 
optimization  process  has  the  added  benefit  of  being  able  to 
easily  resolve  a  location  ambiguity  that  exists  within  the 
system.  Last,  experimental  results  are  presented  using  data 
from  the  2009  data  collection.  Ultimately,  it  is  shown  that  the 
SearchLight  algorithm  enables  precision  radar  geolocation 
using  distributed  receivers  without  the  need  for  precise  timing 
synchronization. 


2.  The  SearchLight  Method 

The  SearchLight  concept  makes  use  of  the  fact  that  a  sta¬ 
tionary,  circularly  scanning  radar  transmitter  will  illuminate 
a  receiver  at  a  regular  interval.  If  one  envisions  two  receivers 
distributed  about  a  circularly  scanning  radar,  the  change  in 
time  of  illumination  between  the  receivers  is  equal  to  the  an¬ 
gle  the  transmitter  sweeps  out  between  the  receivers  divided 
by  the  scan  rate  of  the  transmitter.  This  can  be  expressed  by 
the  following  equation. 


t-j  —  ti  =  —  cos 
u> 


-l 


/(ri_r).(r._r)\ 

V  \ri  —  r\  \rj  —  r|  ) 


(1) 


In  practice,  the  position  of  the  transmitter,  given  here  by 
r,  is  the  unknown  quantity  that  must  be  derived  along  with 
the  scan  rate,  w,  given  a  series  of  differenced  central  times 
of  arrival  (DCTOA),  i,:  and  positions  at  these  times,  r,. 
While  the  position  vector  may  be  expressed  in  either  a  local 
vertical  local  horizontal  coordinate  system  or  in  an  Earth 
fixed  coordinate  system,  the  implicit  assumption  is  that  the 
radar  transmitters  and  receivers  do  not  vary  significantly  from 
the  local  tangent  plane.  The  central  time  of  arrival  for  a 
receiver  is  defined  as  the  time  of  arrival  of  maximum  signal 
amplitude  (corresponding  to  the  arrival  of  the  center  of  the 
main  radar  beam). 


The  method  for  determining  the  central  time  of  arrival  is  de¬ 
pendent  on  radar  signal  characteristics  as  well  as  the  number 
of  active  radar  in  the  collection  area.  When  there  are  multiple 
active  radar  being  collected  against,  each  signal  must  be 
deinterleaved.  There  are  several  methods  for  accomplishing 
this,  all  of  which  rely  on  the  specific  characteristics  of  the 
radar  to  properly  sort  the  signals.  Reference  [2]  contains  a 
general  overview  of  available  techniques.  Once  each  radar 
signal  has  been  isolated,  the  central  time  of  arrival  may  be 
calculated.  For  pulsed  radar,  this  time  may  be  estimated  by 
fitting  a  quadratic  curve  to  the  center  of  received  radar  pulses 
and  calculating  the  time  corresponding  peak  of  the  quadratic 
for  each  radar  pass.  A  similar  technique  may  be  applied  to 
continuous  wave  radar. 


It  should  be  noted  that  the  angle  calculated  in  Equation  1  is  by 
definition  of  arc  cosine  between  0  and  tt;  therefore,  a  check  is 
necessary  not  only  to  ensure  the  correct  quadrant  is  resolved, 
but  also  that  the  calculated  angle  between  receivers  is  indeed 
the  angle  the  radar  has  swept  out  between  detection  events. 


Figure  1.  Diagram  Showing  Example  of  Relative  Geometry 
of  Two  Fixed  Receivers  to  One  Fixed  Radar  Transmitter 


Consider  Figure  1.  As  the  radar  transmitter  sweeps  in  a 
counter-clockwise  direction,  the  angle  it  sweeps  between 
consecutive  receiver  illumination  events  will  alternate  be¬ 
tween  9  and  27t  —  9.  That  is,  the  angle  that  the  radar 
transmitter  sweeps  from  receiver  one  to  receiver  two  is  9, 
while  the  angle  the  transmitter  sweeps  from  receiver  two  to 
receiver  one  is  2tt—9.  However,  the  calculated  angle  between 
the  receivers  will  always  remain  6.  A  set  of  geometric 
checks  must  be  implemented  to  differentiate  these  cases  in 
order  to  calculate  the  correct  error  residual  and  to  extend  this 
technique  to  arbitrary  receiver  geometry. 

It  can  be  shown  that  angle  the  radar  sweeps  between  detection 
events  falls  into  one  of  three  cases.  As  shown  above,  the 
first  case  occurs  when  the  radar  has  swept  out  exactly  the 
calculated  angle  9  between  detection  events.  In  the  next 
case,  the  radar  has  swept  out  not  9  but  its  explement  (i.e. 
27t  —  9).  These  cases  are  demonstrated  in  Figure  1.  The 
last  possible  case  occurs  for  one  mobile  receiver.  Between 
detection  events,  the  receiver  may  have  displaced  an  angle  9 
with  respect  to  itself.  In  this  case,  the  radar  sweeps  out  one 
complete  revolution  plus  the  additional  angular  displacement, 
which  is  27t  +  9.  This  occurs  when  the  receiver  moves 
with  the  direction  of  radar  sweep:  when  the  receiver  moves 
against  the  radar  sweep  the  angle  is  27t  —  9.  Assuming 
that  system  uncertainty  is  small  with  respect  to  the  time 
between  detection  events,  the  correct  DCTOA  error  residual 
is  implemented  as  the  following. 


9  =  cos  1 


/  (■ rt-r )  ■  (r.i+1  —  r)\ 
V  \ri-r\  \rl+1  -  r\  ) 


e  =  DCTOAobser?)  ed  DCTOA  computed 
DCTOA  observed 

DCTOAof,ser„ed  -  (27t  -  9)/u 
pCTOhobserved  -  (27t  +  0)/w 


(2) 

(3) 

(4) 


To  reiterate,  for  each  measurement  a  calculated  difference 
central  time  of  arrival  will  be  computed  based  on  the  cur¬ 
rent  estimate  of  the  radar  transmitter’s  position  and  scan 
rate.  This  calculated  DCTOA  must  be  checked  against  the 
observed  DCTOA  to  ensure  that  the  correct  geometric  case 
is  identified.  From  here,  the  error  residual  is  computed  as 
the  value  that  is  the  minimum  of  the  three  geometric  cases 
discussed  above  and  shown  in  Equation  4.  This  method  is 
valid  as  measurement  error  is  at  least  an  order  of  magnitude 
lower  than  time  difference  between  the  three  geometrical 
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cases.  It  should  be  noted  that  without  making  this  correction, 
the  algorithm  is  constrained  to  geometric  cases  where  each 
receiver  is  180°  or  less  apart  from  each  other  with  respect  to 
the  radar  transmitter.  x  = 


(5) 


Note  that  this  relationship  holds  true  regardless  of  receiver 
configuration.  The  benefit  of  this  method  is  that  data  from 
any  number  of  receivers  in  any  distribution  may  be  used  to 
estimate  the  radar  transmitter’s  location  without  additional 
impact  upon  the  SearchLight  algorithm.  Indeed,  the  only 
requirement  is  that  sufficient  change  in  receiver  bearing  with 
respect  to  the  transmitter  is  observed;  thus,  a  multitude  of 
receiver  platforms  and  configurations  are  supported. 

Further,  it  can  be  seen  that  information  regarding  specific 
transmitter  characteristics  is  not  necessary  a  priori.  This 
gives  SearchLight  an  advantage  of  competing  geolocation 
techniques  such  as  geolocation  via  doppler  shift  observation, 
where  the  emitted  signal  frequency  must  be  known.  Also, 
only  the  arrival  of  the  center  of  the  main  radar  beam  must  be 
isolated  when  multiple  transmitters  are  present,  rather  than 
the  arrival  of  each  pulse.  This  may  decrease  the  complexity 
in  the  signal  deinterleaving  process  and  allow  the  SearchLight 
algorithm  to  be  used  in  applications  where  specific  signal 
characteristics  are  not  known  a  priori. 


1  rrr"1  (  ^  ~  ^  '  ^i+1  ~  r)  \ 
w  V  \ri  ~  r |  \ri+i  —  r\  ) 

(6) 

h(x,ri,ri+ 1) 

(7) 

To  linearize  the  measurement  equation,  the  first  order  vari¬ 
ation  of  its  Taylor  series  expansion  about  a  known  constant 
state  x()Can  be  considered.  So  then,  the  following  equations 
hold. 


y  =  h{x ,  n,  ri+ 1)  -  h(x0,  n,  ri+1) 
y  ss  H  [x  -  jc0] 
dh(x,ri,ri+1) 

H-  a* 


(8) 

(9) 

(10) 


Next,  SearchLight  may  be  used  to  geolocate  Continuous 
Wave  (CW)  radar  transmitters.  Many  traditional  techniques, 
such  as  TDOA,  use  the  arrival  of  pulses  to  define  a  measure¬ 
ment  event  and  cannot  be  easily  implemented  to  geolocate 
CW  radar  transmitters. 

Last,  it  can  be  seen  from  Equation  1  that  this  observation 
is  much  less  sensitive  to  clock  error  than  traditional  time 
difference  of  arrival  (TDOA)  methods.  In  TDOA  methods, 
the  resulting  position  error  due  to  clock  error  is  proportional 
to  the  speed  or  light;  however,  for  this  DCTOA  method, 
the  position  error  is  only  proportional  to  the  scan  rate  of 
the  radar  transmitter,  which  is  several  orders  of  magnitude 
slower.  Thus,  much  larger  clock  errors  are  tolerable.  This 
fact  enables  this  method  to  be  used  with  distributed  receivers 
without  the  need  for  precise  timing  calibration.  Additionally, 
the  requirements  for  data  transfer  with  this  algorithm  is 
much  lower  than  traditional  methods.  Traditional  techniques 
require  data  measurements  to  be  taken  at  the  pulse  rate  of 
the  radar,  while  this  technique  only  requires  the  sharing  of 
data  taken  at  the  scan  rate  of  the  radar  transmitter.  This 
typically  results  in  the  need  for  one  order  of  magnitude  less 
data  transfer  between  receivers  to  formulate  a  geolocation 
solution. 


In  the  weighted  least  squares  (WLS)  formulation  the  goal  is 
to  estimate  the  state,  x°.  that  minimizes  the  following  cost 
functional,  which  is  the  sum  of  the  squares  of  the  error. 

J(x)  =  l/2eTWe  (11) 

=  1/2  (y  -  Hxf  W  (y  -  Hx)  (12) 

Here,  IE  is  a  positive  definite  weighting  matrix  which  is  gen¬ 
erally  set  to  the  inverse  of  the  measurement  noise  covariance. 
This  has  the  well-known  solution  shown  below. 

a;0  =  ( HtWH ) _1  HTWy  (13) 


This  quantity  is  often  described  in  terms  of  the  covariance 
matrix,  P. 

P=(HTWHy 1  (14) 

x°  =  PHTWy  (15) 


3.  Solving  the  Linearized  System 

Given  multiple  DCTOA  measurements,  a  stationary,  circu¬ 
larly  scanning  radar  transmitter’s  position  and  scan  rate  may 
be  estimated.  It  is  clear  that  the  measurement  equation  (Equa¬ 
tion  1)  is  nonlinear,  so  a  linearized  least  squares  solution 
may  be  used.  The  least  squares  criteria  was  first  posed 
by  Gauss.  For  a  more  in  depth  treatment  see  References 
[11]  or  [12].  Here  the  state  that  is  solved  for  is  the  radar 
transmitter’s  position,  r,  and  scan  rate,  u>.  The  measurement 
is  the  DCTOA,  which  is  related  to  the  state  by  Equation  1, 
which  is  repeated  for  completeness. 


As  this  is  a  nonlinear  WLS  problem,  an  initial  state 
estimate, Xq,  must  be  generated  and  the  solution  must  be 
iterated  upon  until  convergence.  An  initial  guess  for  the  radar 
transmitter’s  position  can  be  generated  by  a  rough  grid  search. 
To  generate  an  initial  guess  on  the  radar’s  scan  rate,  the  times 
between  detection  events  for  any  individual  receiver  can  be 
averaged.  This  may  then  be  averaged  between  all  receivers. 
This  time  divided  by  27t  yields  a  fairly  accurate  guess  of  scan 
rate.  This  is  because  the  time  between  detection  events  for 
each  individual  receiver  roughly  corresponds  to  one  complete 
revolution  of  the  radar  beam. 

While  an  analytic  form  of  the  matrix  H  may  be  found 
(this  is  shown  in  the  Appendix),  analysis  has  shown  that 
calculating  numerical  partial  derivatives  generally  yield  better 
results.  Reference  [8]  contains  an  alternate  formulation  of 
this  solution. 
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4.  Solving  the  Nonlinear  System 

As  the  system  is  nonlinear,  it  would  be  beneficial  to  estimate 
the  radar  transmitter’s  position  and  scan  rate  in  a  nonlinear 
fashion.  Additionally,  there  is  an  issue  that  makes  nonlinear 
analysis  preferable.  It  can  be  seen  from  Equation  1  that 
any  given  DCTOA  can  be  generated  by  estimating  the  radar 
transmitter  at  the  true  position  and  true  scan  direction  or  at 
a  position  mirrored  about  the  line  joining  the  measurement 
locations  and  having  the  opposite  scan  direction.  This  false 
location  will  be  referred  to  as  a  “ghost  point.”  This  ghost 
point  will  move  as  the  relative  geometry  of  the  receivers 
change;  however,  the  area  in  the  vicinity  of  the  ghost  point 
will  represent  a  local  minimum  of  the  sum  of  the  squared 
error.  This  may  cause  the  WLS  to  erroneously  converge. 
Thus,  a  method  for  locating  a  global  minimum  (despite  the 
existence  of  local  minimums)  is  necessary.  An  example  of 
this  phenomenon  can  be  seen  in  Figure  2  below.  As  can  be 
seen,  there  is  a  distinct  local  minimum  corresponding  to  the 
ghost  point  mirrored  about  the  mobile  receiver  position. 


accelerated  along  three  directions:  along  its  current  direction, 
towards  its  best  known  individual  solution  and  towards  the 
best  known  global  solution.  This  process  is  iterated  until  a 
stopping  condition  is  met,  which  is  typically  a  fixed  number 
of  iterations. 

The  equations  in  this  process  are  defined  in  the  following 
paragraphs.  Here,  the  objective  function  is  defined  as  the 
sum  of  the  squares  of  the  residual  error:  this  is  a  function 
of  transmitter  position  and  transmitter  scan  rate. 


N 

J(x)  =  55 v 1  ~  hl(x,ri-i .n))2  (16) 

i—2 

Here,  yl  and  h(x.  r  ,  .  r?+1 )’  are  the  ith  observed  DCTOA  and 
calculated  (via  Equation  1  and  Equation  4)  DCTOA,  which 
are  summed  over  M  observations. 


Figure  2.  Plot  Showing  Error  Contours  of  the  Cost  Function 
with  Relative  Receiver  Geometery  for  FI  Dataset 

In  order  to  meet  these  goals,  a  particle  swarm  optimizer 
(PSO)  was  implemented.  PSO  was  first  theorized  in  Ref¬ 
erences  [13]  and  [14].  For  a  more  in  depth  discussion  see 
Reference  [15],  while  Reference  [16]  provides  an  excellent 
overview  of  the  process.  The  problem  must  be  posed  as 
finding  a  minimal  value  of  a  cost  function  given  a  set  of  input 
parameters.  This  applies  exactly  to  the  SearchLight  problem: 
the  sum  of  the  squares  of  the  residuals  may  be  minimized  by 
finding  the  correct  transmitter  location  and  scan  rate.  From 
here,  the  search  space  for  the  solution  must  be  defined.  In 
terms  of  this  application,  this  means  that  appropriate  bounds 
need  to  be  placed  on  the  location  as  well  as  the  scan  rate  for 
the  receiver.  In  this  execution  of  the  PSO,  the  initial  bounds 
on  the  radar  transmitter’s  position  have  been  as  large  as  two 
degrees  in  latitude  and  longitude.  The  initial  bounds  on  the 
scan  rate  can  either  be  estimated  by  the  end  user  or  taken 
as  within  a  reasonable  percentage  of  the  estimated  scan  rate 
using  the  method  described  above.  Next,  the  search  space 
is  randomly  populated  with  a  number  of  state  guesses  (in 
this  case  receiver  position  and  scan  rate),  known  as  particles. 
These  particles  are  also  given  an  arbitrary  initial  velocity. 
At  each  time  step  the  cost  function  is  evaluated  for  each 
particle,  and  the  a  number  of  best  states  are  tracked.  The  first 
tracked  state  is  known  as  the  individual  best  solution  which 
corresponds  to  state  with  the  lowest  cost  at  any  time  step  on  a 
strictly  per  particle  basis.  So,  for  N  particles,  there  will  be  N 
individual  best  states  tracked.  The  next  tracked  state  is  known 
as  the  global  best  solution  and  is  the  state  with  the  lowest  cost 
for  any  particle  at  any  time  step.  There  is  typically  only  one 
global  best  solution.  From  here,  each  particle  is  stochastically 


Now,  an  initial  population  of  particles  can  be  created.  This 
is  done  by  uniformly,  randomly  populating  the  search  space 
with  N  number  of  initial  states.  This  is  accomplished  by 
invoking  the  following  equation  for  the  jth  particle. 


mm  d"  r(0,  1)  d 

d  —  X  max 


(17) 

(18) 


Xj  is  the  state  of  the  jth  particle,  Xmax  is  the  upper  bound 
of  the  search  space,  Xrnin  is  the  lower  bound  for  the  search 
space  and  r(0, 1)  signifies  a  random  number  with  uniform 
distribution  between  0  and  1. 

The  following  process  is  then  iterated  until  a  suitable  stopping 
condition  is  reached.  First,  evaluate  the  objective  function  for 
each  particle. 

Yj  =  J(Xj)  (19) 

Next,  check  to  see  if  this  current  state  represents  either  an 
individual  best  solution  for  the  jth  particle,  denoted  by  -i />  ■, 
or  a  global  best  solution,  represented  by  G.  Here,  only  one 
global  best  solution  will  be  tracked.  If  these  conditions  are 
satisfied  then  the  following  are  implemented. 

*l’j=Xj  if  Yj  <  J(ipj)  (20) 

G  =  Xj  if  Tj  <  J(G)  (21) 

For  M  particles,  the  algorithm  tracks  M  individual  particle 
best  solutions  and  typically  one  global  best  solution.  The  in¬ 
dividual  and  global  best  solutions  are  not  necessarily  updated 
at  every  iteration;  rather,  they  are  only  updated  when  a  better 
individual  particle  solution  or  global  solution  is  found.  This  is 
known  as  a  global  best  topology.  A  formulation  can  be  made 
where  each  particle  is  only  aware  of  the  best  solution  within 
its  group  of  neighbors  as  defined  in  state  space,  as  opposed 
to  each  particle  having  knowledge  of  best  solution  globally. 
This  formulation  is  generally  regarded  as  allowing  the  PSO 
to  be  more  resistant  to  premature  convergence;  however,  it 
requires  more  computation  and  was  not  deemed  necessary  for 
this  work.  Further  discussion  is  available  in  Reference  [17], 

Now  that  all  particles  have  been  evaluated  and  checked  for 
the  individual  or  global  best  solutions,  the  particles  must 
be  moved  to  a  new  position  for  the  next  iteration.  This  is 
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accomplished  in  the  following  manner.  First,  the  velocity  of 
each  particle  is  set  in  the  following  manner. 


v*  =  ^  vT1 

+  cc  (i/>k-Xkj~1) 

+  cs(G-Xk~1)  (22) 

Here,  j  refers  to  the  jth  particle  and  k  refers  to  the  kth  itera¬ 
tion.  The  initial  particle  velocity  can  be  assigned  arbitrarily. 
In  this  implementation,  the  velocity  direction  was  randomly 
assigned  while  the  magnitude  was  set  to  be  arbitrarily  small. 
As  can  been  seen,  the  particle’s  velocity  is  impacted  by 
three  main  factors:  the  particle’s  previous  velocity,  the  vector 
difference  between  its  individual  best  state  and  its  current 
state,  and  the  vector  difference  between  the  globally  best 
solution  and  its  current  state.  These  three  quantities  are  scaled 
by  inertial,  cj,  cognitive,  cc  and  social,  cs  weights,  which  are 
defined  in  the  following  manner. 


1  +  tt(0,  1) 

C/  =  - 2 - 

cc  =  1.49445  r2  (0,1) 

cs  =  1.49445  r3  (0,1)  (23) 

Notice  each  of  these  three  constants  involve  an  independent 
uniform  random  number,  given  by  rj(0, 1).  Each  of  these 
stochastic  scale  factors  were  chosen  from  suggested  values 
in  the  literature  and  may  be  adjusted  in  a  number  of  ways. 
Both  References  [18]  and  [19]  discuss  the  effects  of  altering 
these  parameters.  Once  the  particle  velocity  is  calculated, 
the  particle  is  propagated  to  its  new  location  in  the  following 
manner. 


X)+l=X)  +  V)  (24) 

This  then  concludes  one  iteration.  This  process  is  then 
repeated  until  an  appropriate  stopping  condition  is  met. 


5.  Confidence  Region  Calculation 

In  many  applications,  it  is  not  sufficient  to  only  estimate  the 
position  of  the  radar  transmitter:  it  is  necessary  to  estimate 
the  size  of  the  confidence  region  about  the  solution.  There 
are  two  methods  that  were  explored. 

First,  a  confidence  ellipse  may  be  calculated  from  the  lin¬ 
earized  system.  Reference  [20]  gives  an  excellent  review  of 
calculating  a  confidence  ellipse  given  a  covariance  matrix. 
The  covariance  matrix,  P,  of  the  linearized  system  is  related 
to  the  confidence  ellipse  through  its  eigenvalues,  A j,  and 
eigenvectors,  vl .  As  the  primary  concern  is  for  the  uncer¬ 
tainty  in  the  radar  transmitter’s  position,  consider  the  reduced 
two  by  two  covariance  matrix,  P2x2  (which  corresponds  to 
the  position  covariance),  and  let  its  eigenvalues  be  Ai  >  A2. 
The  confidence  ellipse  can  be  defined  by  its  semi-major  axis, 
SMA,  its  semi-minor  axis,  SMI,  and  its  orientation  angle. 


4>.  These  are  defined  by  the 

following. 

SMA  = 

\/2A  la2F1p~% 

(25) 

SMI  = 

j2\2a2Fl~a_p 

(26) 

Here,  a2  represents  the  calculated  error  variance,  F  rep¬ 
resents  the  percentage  point  of  the  Fisher  distribution,  a 
represents  the  confidence  level,  p  is  the  number  of  degrees 
of  freedom  and  n  is  the  number  of  data  points.  The  Fisher 
distribution  is  itself  a  ratio  of  Gamma  function  and  generally 
has  numerical  values  available  in  standard  math  tables.  The 
orientation  of  the  ellipse  can  defined  as  the  angle  of  the  first 
eigenvector  (i.e.  the  direction  of  the  semi-major  axis)  with 
respect  to  the  first  coordinate  axis.  This  is  given  by  the 
following. 

</>  =  arctan  (27) 

Here,  vj  represents  the  Ith  component  of  the  first  eigenvector. 
This  method  has  the  advantage  of  being  well  suited  to  a  WLS 
approach  as  the  covariance  is  already  calculated.  However,  it 
makes  the  assumption  that  the  confidence  ellipse  is  centered 
upon  the  solution  position,  which  is  not  necessarily  the  case. 
Additionally,  the  ellipse  is  derived  from  the  covariance  of  the 
linearized  system,  which  may  not  yield  the  correct  result  in  a 
nonlinear  system. 

Thus,  a  nonlinear  method  for  developing  the  confidence  re¬ 
gion  is  needed.  References  [21]  and  [22]  both  discuss  several 
methods  for  developing  a  confidence  region  in  the  context  of 
PSO.  A  general  confidence  region  can  be  expressed  by  the 
formula  for  the  likelihood  method  which  is  as  follows. 

J(x)  -  J(x°)  <  a2pFl~a_p  (28) 

This  equation,  which  was  first  formalized  in  Reference  [23], 
states  the  confidence  region  is  the  locus  of  points  x  in  state- 
space  such  that  the  difference  between  the  cost  function  at 
that  point  and  the  cost  function  at  the  optimal  point  x°  is 
within  a  constant,  which  itself  is  a  function  of  the  error 
variance,  degrees  of  freedom  and  the  associated  percentage 
point  of  the  Fisher  distribution. 

This  definition  of  a  confidence  region  lends  itself  well  to  the 
PSO  process.  A  key  feature  of  PSO  is  the  evaluation  of 
the  cost  functional  at  a  large  sampling  of  points.  Thus,  a 
byproduct  of  this  is  the  ability  to  evaluate  whether  a  point 
satisfies  Equation  28,  and  thereby  create  the  locus  of  points 
that  represents  the  confidence  region. 

This  process  is  implemented  in  several  steps.  First,  within 
the  PSO  algorithm  a  check  can  be  added  for  each  particle 
as  to  whether  it  satisfies  Equation  28.  If  it  is  satisfied, 
the  point  and  value  of  the  cost  function  are  stored  for  post 
processing.  After  the  PSO  has  run  to  completion,  a  check 
must  be  run  to  ensure  each  point  that  has  previously  satisfied 
Equation  28  still  indeed  satisfies  the  equation  (as  the  value 
of  the  optimal  cost  has  no  doubt  been  updated  throughout 
each  iteration).  The  remaining  set  of  points  represents  the 
nonlinear  confidence  region.  A  set  of  points,  however,  is 
difficult  to  express  by  a  few  meaningful  quantities.  Thus, 
it  was  decided  that  the  minimum  volume  containing  ellipse 
(MVCE)  should  be  fit  to  these  points  so  that  the  size  and 
orientation  of  region  may  be  expressed  in  identical  terms  as 
in  the  linearized  confidence  region.  In  order  to  make  this 
process  more  efficient,  the  convex  hull  of  this  set  of  points 
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may  be  found  before  calculating  the  MVCE  (since  the  MVCE 
of  a  set  of  points  is  identical  to  the  MVCE  of  its  convex  hull). 
The  convex  hull  is  defined  as  the  minimum  subset  of  points 
such  that  all  points  are  contained  within  them. 

This  process  has  a  few  benefits  derived  from  the  fact  that  the 
confidence  region  is  explicitly  linked  to  the  evaluation  of  the 
cost  functional.  First,  this  method  is  able  to  identify  when 
a  cost  minimum  represented  by  the  ghost  point  is  within  the 
confidence  region  of  the  true  minimum.  This  can  be  a  quality 
measure  by  the  end  user  to  determine  whether  sufficient  data 
has  been  taken  to  definitively  locate  the  radar  transmitter. 
An  additional  benefit  of  this  method  is  that  the  ellipse  is  not 
assumed  to  be  centered  upon  the  optimal  solution.  This  fact 
can  aid  in  identifying  modeling  problems  by  more  readily 
distinguishing  when  the  true  solution  does  not  lie  within  the 
confidence  region  in  excess  of  what  is  statistically  predicted. 


6.  Testing  Results 

An  exercise  was  conducted  where  a  fixed  receiver  and  a 
moving,  ship-based,  mobile  receiver  were  collecting  data 
while  in  range  of  three  separate  radar  transmitters  to  test 
these  methods.  These  transmitters  had  varied  characteristics. 
These  receivers  collected  a  synchronized  time  tag,  their  own 
position  as  well  as  radar  data.  The  relative  geometries  of 
the  radar  transmitters  and  the  receivers  are  shown  in  Figure 
3  below. 

As  can  be  seen,  the  fixed  radar  transmitters  were  installed 
in  a  coastal  location.  The  ship-board  receiver  traveled  along 
the  coast  so  that  the  effects  of  relative  geometries  could 
be  observed.  Figure  3  also  shows  how  the  data  have  been 
divided  into  subsets,  each  set  ranging  from  10  to  30  minutes 
of  collection  time. 

Both  the  WLS  method  and  the  PSO  algorithm  were  imple¬ 
mented  on  the  dataset.  The  estimated  radar  position  was 
compared  against  the  known  position  for  each  data  segment. 
The  results  of  this  analysis  are  summarized  in  Table  1  below. 
As  an  additional  point  of  comparison,  data  for  a  line  of  bear¬ 
ing  (LOB)  geolocation  process  was  simulated.  The  technique 
uses  the  angle  of  arrival  (AOA)  of  the  signal  and  essentially 
uses  triangulation  to  determine  the  radar  transmitter’s  loca¬ 
tion.  References  [2]  and  [5]  contain  detailed  discussion  of 
this  method.  This  was  done  so  that  an  appropriate  comparison 
can  be  made,  as  the  LOB  technique  also  does  not  require  pre¬ 
cision  time  synchronization  and  may  be  used  with  distributed 
receivers.  Both  the  WLS  and  LOB  results  are  taken  from 
Reference  [9]  to  serve  as  a  point  of  comparison  against  the 
PSO  results. 

As  can  be  seen  the  SearchLight  algorithm  provides  geolo¬ 
cation  results  with  much  higher  accuracy  than  the  line  of 
bearing  method.  The  positions  calculated  via  the  SearchLight 
algorithm  are  on  average  twenty  times  more  accurate. 

Furthermore,  it  is  clear  that  the  PSO  generally  provides  more 
accurate  results  than  the  WLS  method.  The  average  miss 
distance  is  almost  150  meters  lower  for  the  PSO  method. 
In  addition  to  the  more  accurate  results,  the  PSO  method 
requires  much  less  data  preprocessing  than  the  WLS  method. 
The  WLS  requires  not  only  an  initial  guess  for  the  radar 
transmitter’s  position  and  scan  rate,  but  it  also  requires  more 
data  checking  to  ensure  that  it  will  not  converge  to  the  ghost 
point.  The  PSO  method,  on  the  other  hand,  requires  only  a 


2WLS  and  LOB  results  are  from  [9] 


Table  1.  Data  Analysis  Results  Summary  of  Miss  Distances 
in  Meters2 


Dataset 

PSO  Method 

WLS  Method 

LOB  Method 

FI 

11.0 

9.3 

926.0 

F2 

38.5 

53.7 

1296.4 

F3 

199.2 

398.2 

6667.2 

F4 

80.4 

196.3 

2592.8 

F5 

160.1 

368.5 

7408.0 

F6 

72.8 

561.2 

14630.8 

F7 

667.0 

822.3 

11112.0 

LI 

72.1 

20.4 

1296.4 

L2 

45.3 

131.5 

9074.8 

N1 

136.5 

785.2 

1852.0 

N2 

732.3 

724.1 

1481.6 

N3 

550.9 

492.6 

2963.2 

Min 

11.0 

9.3 

926.0 

Max 

732.3 

822.3 

14630.8 

Avg 

230.5 

380.3 

5108.4 

rough  search  space  definition  and  is  able  to  uniformly  process 
data  with  little  intervention.  And  since  the  WLS  method 
requires  additional  preprocessing,  the  run  times  for  both 
methods  are  comparable,  even  though  a  PSO  is  traditionally 
more  computationally  intensive. 

Next,  the  confidence  region  of  these  solutions  can  be  calcu¬ 
lated  from  both  the  linearized  and  nonlinear  methods.  These 
regions  can  be  constructed  for  a  given  confidence  level  (e.g. 
95%)  and  represent  where  the  estimated  solution  may  lie  if 
more  or  varied  data  is  used  to  compute  the  solution.  It  should 
be  noted  that  this  does  not  necessarily  mean  that  the  true 
solution  will  lie  within  the  confidence  region  at  that  given 
probability. 

Confidence  ellipses  were  calculated  at  a  99.9%  confidence 
level  in  both  the  linearized  and  nonlinear  method  as  described 
above  for  each  dataset.  It  was  also  noted  whether  the  confi¬ 
dence  region  contained  the  true  solution  or  not.  These  results 
are  presented  in  Table  2  below. 

It  can  be  seen  from  those  results  that  nonlinear  method 
provides  confidence  regions  that  are  much  smaller,  but  only 
contain  the  solution  for  half  of  the  cases  while  the  linearized 
solution  provides  confidence  regions  that  are  much  larger  but 
contain  the  true  solution  for  all  datasets.  At  this  chosen 
confidence  level,  it  would  be  expected  that  each  confidence 
region  would  approximately  always  contain  the  true  solution 
if  the  system  is  correctly  modeled.  This  indicates  that  system 
modeling  can  be  improved.  Indeed,  upon  examination  of  a 
contour  plot  of  the  error  residuals,  it  can  be  seen  that  the 
true  solution  is  not  located  within  a  minimum  error  residual 
location. 

A  plot  of  the  results  for  dataset  FI  are  shown  below  in  Figure 
4.  As  can  be  seen,  the  estimated  radar  transmitter  location  is 
within  1 1  meters  of  the  true  transmitter  location.  The  99.9% 
nonlinear  confidence  ellipse  is  plotted  as  well  for  a  point  of 
reference.  Figure  5  shows  a  plot  of  the  results  for  the  F3 
dataset.  These  results  show  the  case  where  the  true  solution 
is  not  contained  within  the  confidence  ellipse.  It  should  be 
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Figure  3.  Plot  Showing  Relative  Positions  of  Radar  Transmitters  Versus  Receivers  Over  Collection  Period 


Table  2.  Results  Summary  for  Confidence  Region  Analysis 


Dataset 

Nonlinear  Method 

SMA  [m]  SMI  [m]  Contained 

Linearized  Method 

SMA  [m]  SMI  [m]  Contained 

FI 

332.10 

163.34 

In 

288.81 

139.45 

In 

F2 

93.09 

20.36 

Out 

274.02 

57.37 

In 

F3 

1286.39 

28.64 

Out 

7380.09 

80.43 

In 

F4 

313.16 

22.11 

Out 

1660.67 

80.14 

In 

F5 

570.94 

28.49 

Out 

4470.12 

99.37 

In 

F6 

2635.31 

45.32 

Out 

21014.90 

138.75 

In 

F7 

5259.61 

126.09 

In 

50190.45 

315.56 

In 

LI 

1872.81 

202.33 

In 

3584.14 

367.26 

In 

L2 

1749.42 

116.66 

In 

7644.55 

216.75 

In 

N1 

19217.66 

1988.07 

In 

24307.59 

345.76 

In 

N2 

349.64 

155.37 

Out 

1200.12 

508.43 

In 

N3 

7769.74 

6073.10 

In 

30136.18 

72.13 

In 

noted  that  the  confidence  ellipse  extends  beyond  the  bounds 
of  the  plot  and  is  indeed  elliptical.  It  can  be  seen  that  the  true 
transmitter  location  does  not  correspond  to  a  location  with 
minimal  error,  suggesting  that  modeling  errors  are  present. 


7.  Conclusions 

A  novel  method  for  geolocation  of  a  circularly  scanning  radar 
transmitter  based  on  observing  times  between  detection  of 
a  distributed  set  of  receivers  was  introduced.  This  method 
has  the  benefits  of  minimum  complexity  and  minimal  timing 
requirements.  This,  in  turn,  enables  the  technology  to  be  used 
on  a  variety  of  platforms  in  a  multitude  of  configurations 
without  the  need  for  precise  timing  calibration  or  high  data 
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36.932 


Figure  4.  Plot  Showing  Nonlinear  Analysis  Results  for 
Dataset  FI  with  Error  Contours 


Longitude 

Figure  5.  Plot  Showing  Nonlinear  Analysis  Results  for 
Dataset  F3  with  Error  Contours 


transfer  rates.  Both  linear  and  nonlinear  methods  were 
explored  for  solving  this  system.  Using  a  particle  swarm 
optimizer  enabled  precise  geolocation,  while  being  tolerant 
of  the  nonlinearities  of  the  solution  space.  An  analysis  of 
the  confidence  region  associated  with  the  solution  showed 
that  the  uncertainty  in  the  position  of  solution  for  the  radar 
transmitter  is  relatively  small;  however,  it  also  suggests  that 
modeling  improvements  may  be  made. 

With  this  in  mind,  future  work  includes  incorporating  a 
scan  rate  drift  term  in  the  estimation  process.  This  requires 
relatively  low  effort  for  implementation  in  the  PSO  and  will 
allow  a  more  accurate  result  given  a  non-constant  scan  rate. 
In  addition,  incorporating  angle  of  arrival  (AOA)  information 
within  the  estimation  process  should  increase  the  accuracy 
of  the  results  while  maintaining  minimal  complexity.  Uncer¬ 
tainty  in  AOA  measurements  are  generally  orthogonal  to  the 
uncertainty  in  the  DCTOA  measurement.  This  should  allow 
for  increased  accuracy. 


domain  awareness,  are  dependent  upon  the  use  of  a  network 
of  disparate  sensors.  This  method  will  enable  these  sensors 
to  perform  precision  radar  geolocation  without  the  need  for 
expensive  precision  timing  calibration  and  high  data  transfer 
rates. 


Appendix 

The  partial  derivatives  of  the  nonlinear  measurement  given 
in  Equation  1  may  be  calculated  as  follows.  Consider  the 
following  system. 


r 

u ) 

1  cos_,  t (r,-r).(r^r)\ 
w  V  lri  ~  r\  lr2  -  r\  ) 


(29) 

(30) 

(31) 


Now  assume  that  positions  are  restricted  to  the  local  tangent 
plane  and  let  pi  =  rt  —  r.  So  then,  the  following  are  true. 


x  = 


y  =  —  cos 
u> 


-l 


((Pl)‘(P2)\ 

V  IP1IIP2I  ) 


(32) 

(33) 

(34) 


So  then  the  partial  derivatives  of  h  with  respect  to  the  first 
position  state,  x  is  as  follows. 


dh  1  du 

dx  ^  oj  -  ^2  dx 
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du 

dx 


Pi  '  P2 
\Pi\ |p2| 

-Pi(l)  -P2(l) 

IpiI  Ip2| 


+  Pl(!) 


Pi  •  p2 

IpiI3  lp2l 


+  p2(l) 


Pi  •  p2 
lplllp2l3 


(35) 

(36) 

(37) 

(38) 

(39) 


Here,  p(i)  represents  the  ith  component  of  p  vector.  The 
partial  derivative  of  h  with  respect  to  the  second  position 
state,  y,  is  as  follows 


dh  _  1  du 

dy  Twy/1  -72  dy 


du  _  -pi{2)  -  p2{ 2) 

dy  IpiIIp2| 


+  Pi (2) 


Pi  '  P2 

IpiI3  Ip2I 


+  p2(2) 


Pi  '  P2 

|Pil|P2|3 


(40) 

(41) 

(42) 

(43) 


Finally,  the  partial  derivative  of  h  with  respect  to  the  scan 
rate,  w,  is  as  follows. 


Ultimately,  the  technology  shows  enormous  potential  in  a 
variety  of  fields.  The  future  of  some  fields,  such  as  maritime 


dh  1 
—  =  T — 9  arccos 
du; 


Pi  ■  P2  \ 

|piiip2iy 


(44) 


The  sign  ambiguity  on  these  terms  is  related  to  the  ambiguity 
presented  in  Equation  4.  The  partial  derivatives  are  negative 
for  the  case  where  the  radar  has  swept  through  6*  or  27 r  +  6, 
while  the  partial  derivatives  are  positive  for  the  case  where 
the  radar  has  swept  out  2n  —  0. 
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